Housing Prices
Ozan Aygun
11/26/2017
- 1 Introduction and summary
- 2 Loading the data
- 3 Summarizing the data
- 4 Developing expectations from data
- 4.1 Continuous features: watching for strong predictors and potential collinearities
- 4.2 Categorical predictors: feature engineering to create dummy variables
- 4.3 Missing value imputation
- 4.4 GarageYrBlt
- 4.5 LotFrontage
- 4.6 MasVnrArea: Masonry veneer area in square feet
- 4.7 Missing data remaining in the test set
- 5 Predictive modeling using the training set
1 Introduction and summary
2 Loading the data
train <- read.csv("train.csv"); test <- read.csv("test.csv")3 Summarizing the data
summary(train)## Id MSSubClass MSZoning LotFrontage
## Min. : 1.0 Min. : 20.0 C (all): 10 Min. : 21.00
## 1st Qu.: 365.8 1st Qu.: 20.0 FV : 65 1st Qu.: 59.00
## Median : 730.5 Median : 50.0 RH : 16 Median : 69.00
## Mean : 730.5 Mean : 56.9 RL :1151 Mean : 70.05
## 3rd Qu.:1095.2 3rd Qu.: 70.0 RM : 218 3rd Qu.: 80.00
## Max. :1460.0 Max. :190.0 Max. :313.00
## NA's :259
## LotArea Street Alley LotShape LandContour
## Min. : 1300 Grvl: 6 Grvl: 50 IR1:484 Bnk: 63
## 1st Qu.: 7554 Pave:1454 Pave: 41 IR2: 41 HLS: 50
## Median : 9478 NA's:1369 IR3: 10 Low: 36
## Mean : 10517 Reg:925 Lvl:1311
## 3rd Qu.: 11602
## Max. :215245
##
## Utilities LotConfig LandSlope Neighborhood Condition1
## AllPub:1459 Corner : 263 Gtl:1382 NAmes :225 Norm :1260
## NoSeWa: 1 CulDSac: 94 Mod: 65 CollgCr:150 Feedr : 81
## FR2 : 47 Sev: 13 OldTown:113 Artery : 48
## FR3 : 4 Edwards:100 RRAn : 26
## Inside :1052 Somerst: 86 PosN : 19
## Gilbert: 79 RRAe : 11
## (Other):707 (Other): 15
## Condition2 BldgType HouseStyle OverallQual
## Norm :1445 1Fam :1220 1Story :726 Min. : 1.000
## Feedr : 6 2fmCon: 31 2Story :445 1st Qu.: 5.000
## Artery : 2 Duplex: 52 1.5Fin :154 Median : 6.000
## PosN : 2 Twnhs : 43 SLvl : 65 Mean : 6.099
## RRNn : 2 TwnhsE: 114 SFoyer : 37 3rd Qu.: 7.000
## PosA : 1 1.5Unf : 14 Max. :10.000
## (Other): 2 (Other): 19
## OverallCond YearBuilt YearRemodAdd RoofStyle
## Min. :1.000 Min. :1872 Min. :1950 Flat : 13
## 1st Qu.:5.000 1st Qu.:1954 1st Qu.:1967 Gable :1141
## Median :5.000 Median :1973 Median :1994 Gambrel: 11
## Mean :5.575 Mean :1971 Mean :1985 Hip : 286
## 3rd Qu.:6.000 3rd Qu.:2000 3rd Qu.:2004 Mansard: 7
## Max. :9.000 Max. :2010 Max. :2010 Shed : 2
##
## RoofMatl Exterior1st Exterior2nd MasVnrType MasVnrArea
## CompShg:1434 VinylSd:515 VinylSd:504 BrkCmn : 15 Min. : 0.0
## Tar&Grv: 11 HdBoard:222 MetalSd:214 BrkFace:445 1st Qu.: 0.0
## WdShngl: 6 MetalSd:220 HdBoard:207 None :864 Median : 0.0
## WdShake: 5 Wd Sdng:206 Wd Sdng:197 Stone :128 Mean : 103.7
## ClyTile: 1 Plywood:108 Plywood:142 NA's : 8 3rd Qu.: 166.0
## Membran: 1 CemntBd: 61 CmentBd: 60 Max. :1600.0
## (Other): 2 (Other):128 (Other):136 NA's :8
## ExterQual ExterCond Foundation BsmtQual BsmtCond BsmtExposure
## Ex: 52 Ex: 3 BrkTil:146 Ex :121 Fa : 45 Av :221
## Fa: 14 Fa: 28 CBlock:634 Fa : 35 Gd : 65 Gd :134
## Gd:488 Gd: 146 PConc :647 Gd :618 Po : 2 Mn :114
## TA:906 Po: 1 Slab : 24 TA :649 TA :1311 No :953
## TA:1282 Stone : 6 NA's: 37 NA's: 37 NA's: 38
## Wood : 3
##
## BsmtFinType1 BsmtFinSF1 BsmtFinType2 BsmtFinSF2
## ALQ :220 Min. : 0.0 ALQ : 19 Min. : 0.00
## BLQ :148 1st Qu.: 0.0 BLQ : 33 1st Qu.: 0.00
## GLQ :418 Median : 383.5 GLQ : 14 Median : 0.00
## LwQ : 74 Mean : 443.6 LwQ : 46 Mean : 46.55
## Rec :133 3rd Qu.: 712.2 Rec : 54 3rd Qu.: 0.00
## Unf :430 Max. :5644.0 Unf :1256 Max. :1474.00
## NA's: 37 NA's: 38
## BsmtUnfSF TotalBsmtSF Heating HeatingQC CentralAir
## Min. : 0.0 Min. : 0.0 Floor: 1 Ex:741 N: 95
## 1st Qu.: 223.0 1st Qu.: 795.8 GasA :1428 Fa: 49 Y:1365
## Median : 477.5 Median : 991.5 GasW : 18 Gd:241
## Mean : 567.2 Mean :1057.4 Grav : 7 Po: 1
## 3rd Qu.: 808.0 3rd Qu.:1298.2 OthW : 2 TA:428
## Max. :2336.0 Max. :6110.0 Wall : 4
##
## Electrical X1stFlrSF X2ndFlrSF LowQualFinSF
## FuseA: 94 Min. : 334 Min. : 0 Min. : 0.000
## FuseF: 27 1st Qu.: 882 1st Qu.: 0 1st Qu.: 0.000
## FuseP: 3 Median :1087 Median : 0 Median : 0.000
## Mix : 1 Mean :1163 Mean : 347 Mean : 5.845
## SBrkr:1334 3rd Qu.:1391 3rd Qu.: 728 3rd Qu.: 0.000
## NA's : 1 Max. :4692 Max. :2065 Max. :572.000
##
## GrLivArea BsmtFullBath BsmtHalfBath FullBath
## Min. : 334 Min. :0.0000 Min. :0.00000 Min. :0.000
## 1st Qu.:1130 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:1.000
## Median :1464 Median :0.0000 Median :0.00000 Median :2.000
## Mean :1515 Mean :0.4253 Mean :0.05753 Mean :1.565
## 3rd Qu.:1777 3rd Qu.:1.0000 3rd Qu.:0.00000 3rd Qu.:2.000
## Max. :5642 Max. :3.0000 Max. :2.00000 Max. :3.000
##
## HalfBath BedroomAbvGr KitchenAbvGr KitchenQual
## Min. :0.0000 Min. :0.000 Min. :0.000 Ex:100
## 1st Qu.:0.0000 1st Qu.:2.000 1st Qu.:1.000 Fa: 39
## Median :0.0000 Median :3.000 Median :1.000 Gd:586
## Mean :0.3829 Mean :2.866 Mean :1.047 TA:735
## 3rd Qu.:1.0000 3rd Qu.:3.000 3rd Qu.:1.000
## Max. :2.0000 Max. :8.000 Max. :3.000
##
## TotRmsAbvGrd Functional Fireplaces FireplaceQu GarageType
## Min. : 2.000 Maj1: 14 Min. :0.000 Ex : 24 2Types : 6
## 1st Qu.: 5.000 Maj2: 5 1st Qu.:0.000 Fa : 33 Attchd :870
## Median : 6.000 Min1: 31 Median :1.000 Gd :380 Basment: 19
## Mean : 6.518 Min2: 34 Mean :0.613 Po : 20 BuiltIn: 88
## 3rd Qu.: 7.000 Mod : 15 3rd Qu.:1.000 TA :313 CarPort: 9
## Max. :14.000 Sev : 1 Max. :3.000 NA's:690 Detchd :387
## Typ :1360 NA's : 81
## GarageYrBlt GarageFinish GarageCars GarageArea GarageQual
## Min. :1900 Fin :352 Min. :0.000 Min. : 0.0 Ex : 3
## 1st Qu.:1961 RFn :422 1st Qu.:1.000 1st Qu.: 334.5 Fa : 48
## Median :1980 Unf :605 Median :2.000 Median : 480.0 Gd : 14
## Mean :1979 NA's: 81 Mean :1.767 Mean : 473.0 Po : 3
## 3rd Qu.:2002 3rd Qu.:2.000 3rd Qu.: 576.0 TA :1311
## Max. :2010 Max. :4.000 Max. :1418.0 NA's: 81
## NA's :81
## GarageCond PavedDrive WoodDeckSF OpenPorchSF EnclosedPorch
## Ex : 2 N: 90 Min. : 0.00 Min. : 0.00 Min. : 0.00
## Fa : 35 P: 30 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.00
## Gd : 9 Y:1340 Median : 0.00 Median : 25.00 Median : 0.00
## Po : 7 Mean : 94.24 Mean : 46.66 Mean : 21.95
## TA :1326 3rd Qu.:168.00 3rd Qu.: 68.00 3rd Qu.: 0.00
## NA's: 81 Max. :857.00 Max. :547.00 Max. :552.00
##
## X3SsnPorch ScreenPorch PoolArea PoolQC
## Min. : 0.00 Min. : 0.00 Min. : 0.000 Ex : 2
## 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.000 Fa : 2
## Median : 0.00 Median : 0.00 Median : 0.000 Gd : 3
## Mean : 3.41 Mean : 15.06 Mean : 2.759 NA's:1453
## 3rd Qu.: 0.00 3rd Qu.: 0.00 3rd Qu.: 0.000
## Max. :508.00 Max. :480.00 Max. :738.000
##
## Fence MiscFeature MiscVal MoSold
## GdPrv: 59 Gar2: 2 Min. : 0.00 Min. : 1.000
## GdWo : 54 Othr: 2 1st Qu.: 0.00 1st Qu.: 5.000
## MnPrv: 157 Shed: 49 Median : 0.00 Median : 6.000
## MnWw : 11 TenC: 1 Mean : 43.49 Mean : 6.322
## NA's :1179 NA's:1406 3rd Qu.: 0.00 3rd Qu.: 8.000
## Max. :15500.00 Max. :12.000
##
## YrSold SaleType SaleCondition SalePrice
## Min. :2006 WD :1267 Abnorml: 101 Min. : 34900
## 1st Qu.:2007 New : 122 AdjLand: 4 1st Qu.:129975
## Median :2008 COD : 43 Alloca : 12 Median :163000
## Mean :2008 ConLD : 9 Family : 20 Mean :180921
## 3rd Qu.:2009 ConLI : 5 Normal :1198 3rd Qu.:214000
## Max. :2010 ConLw : 5 Partial: 125 Max. :755000
## (Other): 9The outcome variable SalePrice is a complete variable, but there are a number of predictors with many missing values. Median house price is $163,000.
par(mfrow = c(1,2))
boxplot(train$SalePrice,main = "SalePrice")
boxplot(log10(train$SalePrice), main = "Log10(SalePrice)")
3.1 Understanding the impact of missing data
Which features have more than 30% missing data?
na.status <- is.na(train)
na.sum <- apply(na.status,2,sum)
names(na.sum) <- colnames(train)
mostly_missing <- which(na.sum > (0.3 * nrow(train)))
na.sum[mostly_missing]## Alley FireplaceQu PoolQC Fence MiscFeature
## 1369 690 1453 1179 14065 variables have substantial amount of missing data. Do they have impact on the house prices?
3.1.1 Alley
Based on the data description: Type of alley access to property
Grvl Gravel
Pave Paved
NA No alley access
library(ggplot2)
ggplot(data = train, aes(x = Alley, y = log(SalePrice), fill= Alley))+
geom_boxplot()
Alley might be important. Gravel homes have significantly lower price. Convert NAs to “NoAlley” to make more sense.
train$Alley <- as.character(train$Alley)
train$Alley[which(is.na(train$Alley))] <- "NoAlley"
train$Alley <- factor(train$Alley)
#Also transform the test set
test$Alley <- as.character(test$Alley)
test$Alley[which(is.na(test$Alley))] <- "NoAlley"
test$Alley <- factor(test$Alley)3.1.2 FireplaceQu
Based on the data description:
FireplaceQu: Fireplace quality
Ex Excellent - Exceptional Masonry Fireplace
Gd Good - Masonry Fireplace in main level
TA Average - Prefabricated Fireplace in main living area or Masonry Fireplace in basement
Fa Fair - Prefabricated Fireplace in basement
Po Poor - Ben Franklin Stove
NA No Fireplaceggplot(data = train, aes(x = FireplaceQu, y = log10(SalePrice), fill= FireplaceQu))+
geom_boxplot()
Let’s look at through a simple linear model:
par(mfrow=c(1,2))
hist(train$SalePrice)
hist(log10(train$SalePrice)) # Looks more gaussian
FirePlaceFit <- lm(log10(train$SalePrice) ~ FireplaceQu, data = train)
summary(FirePlaceFit)##
## Call:
## lm(formula = log10(train$SalePrice) ~ FireplaceQu, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38691 -0.09351 -0.01374 0.09604 0.57966
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.50252 0.03013 182.644 < 2e-16 ***
## FireplaceQuFa -0.28729 0.03959 -7.256 9.82e-13 ***
## FireplaceQuGd -0.17957 0.03106 -5.781 1.08e-08 ***
## FireplaceQuPo -0.40442 0.04469 -9.050 < 2e-16 ***
## FireplaceQuTA -0.21003 0.03126 -6.719 3.58e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1476 on 765 degrees of freedom
## (690 observations deleted due to missingness)
## Multiple R-squared: 0.118, Adjusted R-squared: 0.1134
## F-statistic: 25.59 on 4 and 765 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(FirePlaceFit)
It appears that all levels of this feature also have significantly different mean Sales prices. Again, converting NAs to NoFireplace to be more informative.
train$FireplaceQu <- as.character(train$FireplaceQu)
train$FireplaceQu[which(is.na(train$FireplaceQu))] <- "NoFireplace"
train$FireplaceQu <- factor(train$FireplaceQu)
#Also transform the test set
test$FireplaceQu <- as.character(test$FireplaceQu)
test$FireplaceQu[which(is.na(test$FireplaceQu))] <- "NoFireplace"
test$FireplaceQu <- factor(test$FireplaceQu)3.1.3 PoolQC
Based on the data description:
PoolQC: Pool quality
Ex Excellent
Gd Good
TA Average/Typical
Fa Fair
NA No Pool
library(ggplot2)
ggplot(data = train, aes(x = PoolQC, y = log10(SalePrice), fill = PoolQC))+
geom_boxplot()
Note that actually very few houses have pools. Houses with excellent pool condition have significantly higher sales prices. We will also keep this feature, converting NAs to “NoPool”:
train$PoolQC <- as.character(train$PoolQC)
train$PoolQC[which(is.na(train$PoolQC))] <- "NoPool"
train$PoolQC <- factor(train$PoolQC)
#Also transform the test set
test$PoolQC <- as.character(test$PoolQC)
test$PoolQC[which(is.na(test$PoolQC))] <- "NoPool"
test$PoolQC <- factor(test$PoolQC)3.1.4 Fence
Based on the data description:
Fence: Fence quality
GdPrv Good Privacy
MnPrv Minimum Privacy
GdWo Good Wood
MnWw Minimum Wood/Wire
NA No Fence library(ggplot2)
ggplot(data = train, aes(x = Fence, y = log10(SalePrice), fill = Fence))+
geom_boxplot()
FenceFit <- lm(log10(train$SalePrice) ~ Fence, data = train)
summary(FenceFit)##
## Call:
## lm(formula = log10(train$SalePrice) ~ Fence, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.57500 -0.06747 -0.00520 0.05611 0.72551
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.23634 0.01816 288.272 < 2e-16 ***
## FenceGdWo -0.11851 0.02628 -4.510 9.58e-06 ***
## FenceMnPrv -0.08969 0.02131 -4.210 3.46e-05 ***
## FenceMnWw -0.11305 0.04582 -2.467 0.0142 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1395 on 277 degrees of freedom
## (1179 observations deleted due to missingness)
## Multiple R-squared: 0.08211, Adjusted R-squared: 0.07217
## F-statistic: 8.26 on 3 and 277 DF, p-value: 2.781e-05Interestingly, this feature might have a complex relationship, if any at all, with the response variable. The effects of each of the level deemed as significant, so we will keep this feature at this point as well.
train$Fence <- as.character(train$Fence)
train$Fence[which(is.na(train$Fence))] <- "NoFence"
train$Fence <- factor(train$Fence)
#Also transform the test set
test$Fence <- as.character(test$Fence)
test$Fence[which(is.na(test$Fence))] <- "NoFence"
test$Fence <- factor(test$Fence)3.1.5 MiscFeature
Based on the data description:
MiscFeature: Miscellaneous feature not covered in other categories
Elev Elevator
Gar2 2nd Garage (if not described in garage section)
Othr Other
Shed Shed (over 100 SF)
TenC Tennis Court
NA None
We will deal with the remaining missing value containing features in the next stages of our analysis.
table(train$MiscFeature)##
## Gar2 Othr Shed TenC
## 2 2 49 1ggplot(data = train, aes(x = MiscFeature, y = log10(SalePrice), fill = MiscFeature))+
geom_boxplot()
Again, very few houses have this features, and not sure if these qualities have a real impact on the house value amongst other measured features. Neverthless, I will maintain this feature as well. Just converting NAs to “NoMiscFeature” to be more clear:
train$MiscFeature <- as.character(train$MiscFeature)
train$MiscFeature[which(is.na(train$MiscFeature))] <- "NoMiscFeature"
train$MiscFeature <- factor(train$MiscFeature)
#Also transform the test set
test$MiscFeature <- as.character(test$MiscFeature)
test$MiscFeature[which(is.na(test$MiscFeature))] <- "NoMiscFeature"
test$MiscFeature <- factor(test$MiscFeature)4 Developing expectations from data
Here I will explore the data further, in order to develope a better understanding about the relationships between features and the home prices.
4.1 Continuous features: watching for strong predictors and potential collinearities
Let’s look at the continuous features in the data set:
cont.var <-NULL
for(i in 1:ncol(train)){
if(class(train[,i]) == "integer" | class(train[,i]) == "numeric"){
cont.var <- c(cont.var,i)
}
}4.1.1 Part1: transforming LotArea to make it more gaussian
pairs(log10(train$SalePrice) ~ ., data = train[,cont.var[1:10]], cex = 0.05, pch = 19, col = "navy")
- We need to remove ID, since this is just an identifier for the houses, has no predictive value as expected.
- OverallQual and OverallCond look strong predictors, but might have some collinearity.
- Yearbuilt and YearRemodAdd also look good predictors but most likely contain redundant information.
- For the LOT and SF features, it could be better to look at after log transformation, since the data looks skewed.
par(mfrow=c(3,2))
hist(train$LotArea, breaks = 50, col = "navy");hist(log10(train$LotArea+1),breaks=50,col = "navy")
hist(train$MasVnrArea, breaks = 50,col = "navy");hist(log10(train$MasVnrArea+1),breaks=50,col = "navy")
hist(train$BsmtFinSF1, breaks = 50,col = "navy");hist(log10(train$BsmtFinSF1+1),breaks=50,col = "navy")
Indeed one of these features (LotArea) look more gaussian when it is log10 transformed.
pairs(log10(train$SalePrice) ~ log10(train$LotArea+1) + log10(train$MasVnrArea+1) + log10(train$BsmtFinSF1+1),cex = 0.1, pch = 19, col = "navy")
When transformed, lotArea seems to do a better job in explaining the variability in the SalePrice, other 2 variables however, have substantial 0 values and it is questionable how much this transformation might help.
Therefore, I will only transform the LotArea at this stage and leave the other features as they are, for feature selection steps.
4.1.2 Part2: Area related features generally explain variation in the SalePrice
pairs(log10(train$SalePrice) ~ ., data = train[,cont.var[11:20]], cex = 0.05, pch = 19, col = "navy")
4.1.3 Part3: Size of the garage and number of rooms are good predictors
pairs(log10(train$SalePrice) ~ ., data = train[,cont.var[22:31]], cex = 0.05, pch = 19, col = "navy")
Arguably, these features are all related to the overall area of the house.
4.1.4 Part4: fairly weak predictors
pairs(log10(train$SalePrice) ~ ., data = train[,cont.var[32:38]], cex = 0.05, pch = 19, col = "navy")
Processing the identified features:
# Remove ID variable
train.ID <- train$Id #Hold it in a new variable
train <- train[, -1]
library(dplyr)
train <- dplyr::mutate(train, log10.LotArea = log10(LotArea))
train <- dplyr::select(train, -LotArea)
# Also transform in the test set
test <- dplyr::mutate(test, log10.LotArea = log10(LotArea))
test <- dplyr::select(test, -LotArea)4.2 Categorical predictors: feature engineering to create dummy variables
Now let’s look at the categorical variables to see if certain ones are more powerful than others:
cat.var <-NULL
for(i in 1:ncol(train)){
if(class(train[,i]) == "factor"){
cat.var <- c(cat.var,i)
}
}
length(cat.var)## [1] 43We have 43 categorical variables, some of them have only few levels. This is a challenge in many ways. It is not feasible to look each of them individually, therefore, we will first convert them to binary variables and try to use regularization or decision tree-based approaches to select and use most useful features.
# Write a function that gets a data set and converts all factor variables to dummy variables
convert.to.dummy <- function(data.set){
cat.var <-NULL
temp.data <- data.frame(1:nrow(data.set))
for(i in 1:ncol(data.set)){
if(class(data.set[,i]) == "factor"){
cat.var <- c(cat.var,i)
factor.levels <- levels(data.set[,i]) # Try to find a way to classify NA's as "NO" otherwise they generate problem downstream
# First check if there is any 'NA-level'
if(any(is.na(data.set[,i]))){
dummy.vector = ifelse(is.na(data.set[,i]),1,0)
dummy.vector <- data.frame(dummy.vector)
colnames(dummy.vector)[1] = paste("NO",names((data.set)[i]),sep = ".")
temp.data <- cbind(temp.data,dummy.vector)
}
for(j in seq_along(factor.levels)){ # Then deal with normal factor levels
dummy.vector = ifelse(data.set[,i] == factor.levels[j],1,0)
#Since we already dealt with NAs above
if(any(is.na(dummy.vector))){dummy.vector[is.na(dummy.vector)] <- 0}
dummy.vector <- data.frame(dummy.vector)
colnames(dummy.vector)[1] = paste(names((data.set)[i]),
factor.levels[j],sep = ".")
temp.data <- cbind(temp.data,dummy.vector)
}
}
}
#Remove the original categorical variables from data.set
data.set <- data.set[,-cat.var]
#Add the dummy.variable set
temp.data <- temp.data[,-1] # remove the unnecessary column
data.set <- cbind(data.set,temp.data)
return(data.set)
}
# Keep test.IDs aside
test.ID <- test$Id
# Process the training set
training.processed <- convert.to.dummy(train)
# Process the test set
test.processed <- convert.to.dummy(test)
# Note that not all the levels are present in the test set, therefore we need to make the features same.
training.processed.SalePrice <- training.processed$SalePrice #Keep the original outcome variable
consensus.features1 <- which(colnames(training.processed) %in% colnames(test.processed))
training.processed <- training.processed[,consensus.features1] #Keep the same features as in the test set
consensus.features2 <- which(colnames(test.processed) %in% colnames(training.processed))
test.processed <- test.processed[,consensus.features2] #Keep the same features as in the training set
identical(colnames(test.processed), colnames(training.processed))## [1] TRUE# Same features in both sets!
training.processed <- data.frame(Log.SalePrice = log(training.processed.SalePrice),training.processed) # Add the Log of the response variable
test.processed <- data.frame(ID = test.ID,test.processed) # Add the ID variable4.3 Missing value imputation
We still have missing values in some of our features, which will be a problem down the road.
apply(is.na(training.processed), 2, sum)[apply(is.na(training.processed), 2, sum) != 0]## LotFrontage MasVnrArea GarageYrBlt
## 259 8 81apply(is.na(test.processed), 2, sum)[apply(is.na(test.processed), 2, sum) != 0]## LotFrontage MasVnrArea BsmtFinSF1 BsmtFinSF2 BsmtUnfSF
## 227 15 1 1 1
## TotalBsmtSF BsmtFullBath BsmtHalfBath GarageYrBlt GarageCars
## 1 2 2 78 1
## GarageArea
## 14.4 GarageYrBlt
We have 81 missing cases in our training set, these are the houses actually they don’t have a garage.
plot(y =training.processed$Log.SalePrice, x =training.processed$GarageYrBlt, pch = 19, col = "navy", cex = 0.5)
It looks like this feature has some relationship with the outcome, as also we can generally expect that newer houses are going to be more expensive. Therefore, we should make an attempt to impute the values missing from this feature.
plot(training.processed$YearBuilt, training.processed$GarageYrBlt, pch = 19, col = "navy", cex = 0.5)
As we expected, for most of the houses the year that garage was built is the same as the year the house was built. Therefore, for imputation purpose it is sensible to use the YearBuilt feature to fill these issing values in the GarageYrBuilt feature. While this is not consistent with the fact that these houses do not contain any garage, since we captured that in the categorical features we engineered, at this point imputing the years seems to be more appropriate than removing these data points alltogether:
#Impute in the training set
training.processed$GarageYrBlt[is.na(training.processed$GarageYrBlt)] <- training.processed$YearBuilt[is.na(training.processed$GarageYrBlt)]
#Impute in the test set
test.processed$GarageYrBlt[is.na(test.processed$GarageYrBlt)] <- test.processed$YearBuilt[is.na(test.processed$GarageYrBlt)]4.5 LotFrontage
Let’s look at this feature more closely:
plot(x = training.processed$LotFrontage, y = training.processed$Log.SalePrice, cex = 0.5, col = "navy" , pch = 19)
abline(lm(training.processed$Log.SalePrice ~ training.processed$LotFrontage))
Again this feature seems to have some degree of relationship with the response. We suspect that this feature has relationship with other features that relate to the Lot or the general area of the house:
cor(training.processed$LotFrontage,training.processed$log10.LotArea, use = "complete.obs")## [1] 0.6540587cor(training.processed$LotFrontage,training.processed$GarageArea, use = "complete.obs")## [1] 0.3449967cor(training.processed$LotFrontage,training.processed$MasVnrArea, use = "complete.obs")## [1] 0.1934581# Note that correlation gets better when log10 transformed:
cor(log10(training.processed$LotFrontage), training.processed$log10.LotArea, use = "complete.obs")## [1] 0.7455501summary(lm(training.processed$Log.SalePrice ~ log10(training.processed$LotFrontage) + training.processed$log10.LotArea ))##
## Call:
## lm(formula = training.processed$Log.SalePrice ~ log10(training.processed$LotFrontage) +
## training.processed$log10.LotArea)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.51402 -0.23961 -0.04072 0.25528 1.22485
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.88631 0.20139 44.125 < 2e-16
## log10(training.processed$LotFrontage) 0.30471 0.10638 2.864 0.00425
## training.processed$log10.LotArea 0.65296 0.07572 8.623 < 2e-16
##
## (Intercept) ***
## log10(training.processed$LotFrontage) **
## training.processed$log10.LotArea ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3763 on 1198 degrees of freedom
## (259 observations deleted due to missingness)
## Multiple R-squared: 0.1833, Adjusted R-squared: 0.1819
## F-statistic: 134.4 on 2 and 1198 DF, p-value: < 2.2e-16plot(y = log10(training.processed$LotFrontage), x = training.processed$log10.LotArea, cex = 0.5, col = "navy" , pch = 19)
abline(lm(log10(training.processed$LotFrontage) ~ training.processed$log10.LotArea))
It looks like this feature pretty much collinear with the LotArea feature, especially when log10 transformed. Therefore, we will remove this feature from both training and test data sets.
training.processed <- dplyr::select(training.processed, -LotFrontage)
test.processed <- dplyr::select(test.processed, -LotFrontage)4.6 MasVnrArea: Masonry veneer area in square feet
Looking at more closely to this feature:
plot(x = training.processed$MasVnrArea, y = training.processed$Log.SalePrice,cex = 0.5, col = "navy" , pch = 19)
Note that a large number of these values are zero, which probably reflects the houses with no masonry veener. What happens if we put these aside and look at the rest of the values:
plot(x = log10(training.processed$MasVnrArea[training.processed$MasVnrArea != 0]), y = training.processed$Log.SalePrice[training.processed$MasVnrArea != 0],cex = 0.5, col = "navy" , pch = 19)
There might be ‘some’ relationship with the response, if not too strong. It is also a little confusing since the variable MasVnrType only has these 8 values as “none”, but linked the large amount of MasVnrAre == 0 data to certain MasVnrTypes:
qplot(data = training.processed,x = log10(training.processed$MasVnrArea), y = training.processed$Log.SalePrice, col = train$MasVnrType)
It is more clear now. The “none” group indeed represent the MasVnrArea = 0, NAs however appear as missing values in both cases.
ggplot(data = training.processed,aes(x = train$MasVnrType, y = training.processed$Log.SalePrice, fill = train$MasVnrType))+geom_boxplot()
We note that the SalePrice of the missing cases in the training set is closest to the SalePrice of the houses with Stone Masonry Veener type. It is likely that the area of the veener in some of these houses were not measured, perhaps it is more difficult to do so. This is our assumption.
Based on this assumption, we will merge the NA cases with the Stone group and remove the MasVnrArea feature in both training and test data sets:
training.processed$MasVnrType.None[training.processed$NO.MasVnrType == 1] <- 1
training.processed <- dplyr::select(training.processed, -NO.MasVnrType, -MasVnrArea)
test.processed$MasVnrType.None[test.processed$NO.MasVnrType == 1] <- 1
test.processed <- dplyr::select(test.processed, -NO.MasVnrType, -MasVnrArea)4.7 Missing data remaining in the test set
We still need to deal with few missing data points we have in the test set. Since we need to make SalePrice predictions for each of these cases, we need to assess how we can impute these values at this point:
apply(is.na(test.processed), 2, sum)[apply(is.na(test.processed), 2, sum) != 0]## BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF BsmtFullBath
## 1 1 1 1 2
## BsmtHalfBath GarageCars GarageArea
## 2 1 1Some features regarding Basement and Garage have missing values:
test.processed[which(apply(is.na(test.processed), 1, sum) != 0),which(apply(is.na(test.processed), 2, sum) != 0 )]Fortunately, the missing values are only restricted to 3 house records. It is likely these unusual houses are either lacking basement or garage. Since all these features are numeric, we will simply place zero instead of NA’s for these values. This is of course based on our assumption.
for (i in 1:ncol(test.processed)) {
if(any(is.na(test.processed[,i]))){
test.processed[is.na(test.processed[,i]),i] <- 0
}
}This completes missing value imputation as no missing values left in either training and test data sets.
5 Predictive modeling using the training set
5.1 Using shrinkage methods (regularization)
It would be interesting to perform a linear model selection and see how well it might perform with the new training set we obtained.
5.1.1 Fitting a lasso regression model to data
library(glmnet)
# We create a matrix of predictors
x.training = model.matrix(Log.SalePrice ~ . -1, data = training.processed)
x.testing = model.matrix( ~ . -1, data = test.processed)
y = training.processed$Log.SalePrice
# Fit lasso regression to training data
fit.lasso <- glmnet(x.training,y, family = "gaussian")
plot(fit.lasso, xvar = "lambda", label = TRUE)
plot(fit.lasso, xvar = "dev", label = TRUE)
Lasso is a nice way of reducing the features while still trying to explain the variance and therefore trying to win the bias-variance trade-off. In this case we notice that only 3 features are able to explain 60% of the deviance, while adding 13 features is able to explain 80% of it. It is more feasible to use the restricted model with the few features and their shrunken coefficients.
5.1.2 Choosing the optimal lambda by using cross-validation
It is best to choose the optimal lambda using cross-validation, with the aim of minimizing mean squared error::
#cv.lasso will be the list containing the coefficients and information of our optimal model fit using cross-validation (10 fold by default)
cv.lasso <- cv.glmnet(x.training,y, family = "gaussian")
plot(cv.lasso)
Notice that what we are doing here is actually computationally intense.
- We get very fine grids of lambda values.
- Using each lambda value, we perform 10-fold cross validation:
- We fit a model using the Sum of Squares convex optimization and obtain cross-validated error (average of 10 model fits)
- We continue steps 1 and 2 for all fine grids of lambda values in the range.
The two vertical lines are produced in the plot, the left one marks the model with min error and the middle one shows a little more restricted model 1 standard error away from the minimum error.
When choosing the optimal model, we might prefer to get the more restricted model that is indicated by the second vertical line, which is default in glmnet package, due to the 1 standard error conversion.
The final (optimal) model coefficients we obtained by cross-validation and the non-zero features remaining:
coef(cv.lasso)## 283 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 5.309391e+00
## MSSubClass .
## OverallQual 9.324533e-02
## OverallCond 1.040886e-02
## YearBuilt 1.182768e-03
## YearRemodAdd 1.355597e-03
## BsmtFinSF1 1.629126e-05
## BsmtFinSF2 .
## BsmtUnfSF .
## TotalBsmtSF 5.537271e-05
## X1stFlrSF .
## X2ndFlrSF .
## LowQualFinSF .
## GrLivArea 1.803022e-04
## BsmtFullBath 1.643841e-02
## BsmtHalfBath .
## FullBath .
## HalfBath .
## BedroomAbvGr .
## KitchenAbvGr .
## TotRmsAbvGrd .
## Fireplaces 9.766781e-03
## GarageYrBlt .
## GarageCars 6.852253e-02
## GarageArea .
## WoodDeckSF 1.644715e-05
## OpenPorchSF .
## EnclosedPorch .
## X3SsnPorch .
## ScreenPorch .
## PoolArea .
## MiscVal .
## MoSold .
## YrSold .
## log10.LotArea 1.533081e-01
## MSZoning.C..all. -1.476063e-01
## MSZoning.FV .
## MSZoning.RH .
## MSZoning.RL .
## MSZoning.RM -4.101732e-02
## Street.Grvl .
## Street.Pave .
## Alley.Grvl .
## Alley.NoAlley .
## Alley.Pave .
## LotShape.IR1 .
## LotShape.IR2 .
## LotShape.IR3 .
## LotShape.Reg .
## LandContour.Bnk .
## LandContour.HLS .
## LandContour.Low .
## LandContour.Lvl .
## Utilities.AllPub .
## LotConfig.Corner .
## LotConfig.CulDSac .
## LotConfig.FR2 .
## LotConfig.FR3 .
## LotConfig.Inside .
## LandSlope.Gtl .
## LandSlope.Mod .
## LandSlope.Sev .
## Neighborhood.Blmngtn .
## Neighborhood.Blueste .
## Neighborhood.BrDale .
## Neighborhood.BrkSide .
## Neighborhood.ClearCr .
## Neighborhood.CollgCr .
## Neighborhood.Crawfor 1.844449e-02
## Neighborhood.Edwards .
## Neighborhood.Gilbert .
## Neighborhood.IDOTRR .
## Neighborhood.MeadowV .
## Neighborhood.Mitchel .
## Neighborhood.NAmes .
## Neighborhood.NoRidge .
## Neighborhood.NPkVill .
## Neighborhood.NridgHt 5.162883e-04
## Neighborhood.NWAmes .
## Neighborhood.OldTown .
## Neighborhood.Sawyer .
## Neighborhood.SawyerW .
## Neighborhood.Somerst .
## Neighborhood.StoneBr .
## Neighborhood.SWISU .
## Neighborhood.Timber .
## Neighborhood.Veenker .
## Condition1.Artery .
## Condition1.Feedr .
## Condition1.Norm 4.167761e-03
## Condition1.PosA .
## Condition1.PosN .
## Condition1.RRAe .
## Condition1.RRAn .
## Condition1.RRNe .
## Condition1.RRNn .
## Condition2.Artery .
## Condition2.Feedr .
## Condition2.Norm .
## Condition2.PosA .
## Condition2.PosN -5.019406e-02
## BldgType.1Fam .
## BldgType.2fmCon .
## BldgType.Duplex .
## BldgType.Twnhs .
## BldgType.TwnhsE .
## HouseStyle.1.5Fin .
## HouseStyle.1.5Unf .
## HouseStyle.1Story .
## HouseStyle.2.5Unf .
## HouseStyle.2Story .
## HouseStyle.SFoyer .
## HouseStyle.SLvl .
## RoofStyle.Flat .
## RoofStyle.Gable .
## RoofStyle.Gambrel .
## RoofStyle.Hip .
## RoofStyle.Mansard .
## RoofStyle.Shed .
## RoofMatl.CompShg .
## RoofMatl.Tar.Grv .
## RoofMatl.WdShake .
## RoofMatl.WdShngl .
## Exterior1st.AsbShng .
## Exterior1st.AsphShn .
## Exterior1st.BrkComm .
## Exterior1st.BrkFace .
## Exterior1st.CBlock .
## Exterior1st.CemntBd .
## Exterior1st.HdBoard .
## Exterior1st.MetalSd .
## Exterior1st.Plywood .
## Exterior1st.Stucco .
## Exterior1st.VinylSd .
## Exterior1st.Wd.Sdng .
## Exterior1st.WdShing .
## Exterior2nd.AsbShng .
## Exterior2nd.AsphShn .
## Exterior2nd.Brk.Cmn .
## Exterior2nd.BrkFace .
## Exterior2nd.CBlock .
## Exterior2nd.CmentBd .
## Exterior2nd.HdBoard .
## Exterior2nd.ImStucc .
## Exterior2nd.MetalSd .
## Exterior2nd.Plywood .
## Exterior2nd.Stone .
## Exterior2nd.Stucco .
## Exterior2nd.VinylSd .
## Exterior2nd.Wd.Sdng .
## Exterior2nd.Wd.Shng .
## MasVnrType.BrkCmn .
## MasVnrType.BrkFace .
## MasVnrType.None .
## MasVnrType.Stone .
## ExterQual.Ex .
## ExterQual.Fa .
## ExterQual.Gd .
## ExterQual.TA .
## ExterCond.Ex .
## ExterCond.Fa .
## ExterCond.Gd .
## ExterCond.Po .
## ExterCond.TA .
## Foundation.BrkTil .
## Foundation.CBlock .
## Foundation.PConc .
## Foundation.Slab .
## Foundation.Stone .
## Foundation.Wood .
## NO.BsmtQual .
## BsmtQual.Ex 2.446339e-02
## BsmtQual.Fa .
## BsmtQual.Gd .
## BsmtQual.TA .
## NO.BsmtCond .
## BsmtCond.Fa .
## BsmtCond.Gd .
## BsmtCond.Po .
## BsmtCond.TA .
## NO.BsmtExposure .
## BsmtExposure.Av .
## BsmtExposure.Gd 9.061738e-04
## BsmtExposure.Mn .
## BsmtExposure.No .
## NO.BsmtFinType1 .
## BsmtFinType1.ALQ .
## BsmtFinType1.BLQ .
## BsmtFinType1.GLQ 8.339921e-03
## BsmtFinType1.LwQ .
## BsmtFinType1.Rec .
## BsmtFinType1.Unf -6.834311e-03
## NO.BsmtFinType2 .
## BsmtFinType2.ALQ .
## BsmtFinType2.BLQ .
## BsmtFinType2.GLQ .
## BsmtFinType2.LwQ .
## BsmtFinType2.Rec .
## BsmtFinType2.Unf .
## Heating.GasA .
## Heating.GasW .
## Heating.Grav .
## Heating.Wall .
## HeatingQC.Ex 1.334368e-02
## HeatingQC.Fa .
## HeatingQC.Gd .
## HeatingQC.Po .
## HeatingQC.TA .
## CentralAir.N -6.204003e-02
## CentralAir.Y .
## Electrical.FuseA .
## Electrical.FuseF .
## Electrical.FuseP .
## Electrical.SBrkr .
## KitchenQual.Ex 4.700208e-03
## KitchenQual.Fa .
## KitchenQual.Gd .
## KitchenQual.TA -1.387015e-02
## Functional.Maj1 .
## Functional.Maj2 .
## Functional.Min1 .
## Functional.Min2 .
## Functional.Mod .
## Functional.Sev .
## Functional.Typ .
## FireplaceQu.Ex .
## FireplaceQu.Fa .
## FireplaceQu.Gd .
## FireplaceQu.NoFireplace -2.957736e-02
## FireplaceQu.Po .
## FireplaceQu.TA .
## NO.GarageType .
## GarageType.2Types .
## GarageType.Attchd 6.147442e-03
## GarageType.Basment .
## GarageType.BuiltIn .
## GarageType.CarPort .
## GarageType.Detchd .
## NO.GarageFinish .
## GarageFinish.Fin .
## GarageFinish.RFn .
## GarageFinish.Unf .
## NO.GarageQual .
## GarageQual.Fa .
## GarageQual.Gd .
## GarageQual.Po .
## GarageQual.TA .
## NO.GarageCond .
## GarageCond.Ex .
## GarageCond.Fa .
## GarageCond.Gd .
## GarageCond.Po .
## GarageCond.TA .
## PavedDrive.N -5.958627e-03
## PavedDrive.P .
## PavedDrive.Y .
## PoolQC.Ex .
## PoolQC.Gd .
## PoolQC.NoPool .
## Fence.GdPrv .
## Fence.GdWo .
## Fence.MnPrv .
## Fence.MnWw .
## Fence.NoFence .
## MiscFeature.Gar2 .
## MiscFeature.NoMiscFeature .
## MiscFeature.Othr .
## MiscFeature.Shed .
## SaleType.COD .
## SaleType.Con .
## SaleType.ConLD .
## SaleType.ConLI .
## SaleType.ConLw .
## SaleType.CWD .
## SaleType.New .
## SaleType.Oth .
## SaleType.WD .
## SaleCondition.Abnorml -7.582707e-03
## SaleCondition.AdjLand .
## SaleCondition.Alloca .
## SaleCondition.Family .
## SaleCondition.Normal .
## SaleCondition.Partial .Note that our final model includes ~30 features with non-zero coefficients. We therefore significantly reduced the size of the data set by selecting features through lasso regression. The resulting model explains more than 80% of the deviance in the outcome, which is considerable good for the training set.
We also note that some of the dummy variables remained in the model, while majority of them removed. It seems to be good idea to keep all levels as seperate dummy variables until the regularization, since we don’t know which ones may be important in the final model context.
5.1.3 Making predictions using the best lasso model fit
Let’s first start with the training set:
pred <- predict(cv.lasso, newx = x.training)
#Calculating the root mean squared error
rmse.training = sqrt(mean((training.processed$Log.SalePrice - pred)^2))
rmse.training## [1] 0.1479088#Pearson correlation:
cor.training = cor(training.processed$Log.SalePrice,pred)
plot(x = pred, y = training.processed$Log.SalePrice, cex = 0.5, col = "navy", pch = 19)
text(x = max(pred), y = min(training.processed$Log.SalePrice)+0.3,
labels = paste0("RMSE: ",rmse.training))
Therefore, the RMSE we obtained from the training set is ~15%.
Next we will try to use test set to estimate the house prices using our model.